The Coulomb impurity problem in graphene 
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We address the problem of an unscreened Coulomb charge in graphene, and calculate the local 
density of states and displaced charge as a function of energy and distance from the impurity. 
This is done non-perturbatively in two different ways: (1) solving the problem exactly by studying 
numerically the tight-binding model on the lattice; (2) using the continuum description in terms of 
the 2D Dirac equation. We show that the Dirac equation, when properly regularized, provides a 
qualitative and quantitative low energy description of the problem. The lattice solution shows extra 
features that cannot be described by the Dirac equation, namely bound state formation and strong 
renormalization of the van Hove singularities. 

PACS numbers: 81.05.Uw, 71.55.-i, 25.75.Dw 



Since the isolation of graphene a few years ago [1] 
there is an intense interest in understanding the elec- 
tronic properties of this material. The low energy elec- 
tronic excitations in graphene arc linearly dispersing, 
massless, chiral Dirac fcrmions, described by Dirac cones 
at the edges of the Brillouin zone (BZ) (at the K and K' 
points). Due to the vanishing density of states of these 
Dirac fermions in undoped graphene, the system is very 
sensitive to impurities and defects [2], which control its 
transport properties [3, 4]. This also means that elec- 
trons in graphene screen poorly, and hence the issue of 
unscreened Coulomb interactions becomes paramount in 
the understanding of the experimental data [1] . 

In this paper we contrast the tight-binding approach 
(that we solve exactly with numerical techniques) with 
the continuum approach based on the Dirac equation. 
We show that the latter provides a good qualitative de- 
scription of the problem at low energies, when properly 
regularized. We also show that the Dirac description fails 
at moderate to high energies and at short distances, when 
the lattice description is the only one possible. In this 
case new features, not captured by the Dirac Hamilto- 
nian, emerge. We calculate the local density of states 
(LDOS) and induced charge around a Coulomb impurity 
as a function of energy and distance. These quantities 
are experimentally accessible through scanning tunneling 
spectroscopy (STS). We stress that calculations for im- 
purities with long-range potentials are radically different 
from the ones for short-range forces, which are exactly 
solvable using T-matrix methods [5]. 

Consider the problem of a single Coulomb impurity, 
with charge Ze, placed in the middle of an hexagon of the 
honeycomb lattice [6]. The tight-binding Hamiltonian 
for this problem, with nearest-neighbor hopping only, is 
given by (we use units such that h = 1): 



U = i^(af&j +h.c. 



Ze 2 



E 



b\bi 



(1) 



where m (bi) annihilates an electron at site and sub- 
lattice A (B), t w 2.7 cV is the hopping energy, and 
rf' B is the distance between the carbon atoms and the 



impurity (assumed to be at the origin of the coordinate 
system); £o is the dielectric constant. We calculated nu- 
merically the spectrum of (1) using the methods of exact 
diagonalization and recursion used in ref. [7] for the study 
of short-range unitary scatterers. 

Close to the K point in the BZ we can write an ef- 
fective low energy Hamiltonian for (1) in terms of Dirac 
fermions with a spinor wave function ^(r), whose compo- 
nents represent its weight on each sublattice. The wave 
function obeys the equation: 



up [a ■ p- g/r)^(r) = E *(r), 



(2) 



with vp = iat/2 (« 10 6 m/s) being the Fermi veloc- 
ity, p the 2D momentum operator, <7j the Pauli matrices, 
and g = Ze 2 /(up Eq) is the dimensionless coupling con- 
stant. Henceforth, we shall take a (the C-C distance) 
and up as distance and energy units. Notice that (2) 
does not involve inter-cone scattering because the un- 
screened Coulomb potential is dominated by small mo- 
mentum transfers since, in Fourier space, it behaves like 
1/q and is singular as q — > 0. 

Eq. (2) is separable in cylindrical coordinates. Resort- 
ing to eigenfunctions of the conserved angular momen- 
tum, J z = L z + a z /2 [8], 



(3) 



the radial equation for (2) reads (j = ±1/2, ±3/2, . . .) 



e+g/r -(d r +j/r) 
(dr-j/r) e+g/r 



= M m (r) = 0. (4) 



This equation can be solved by multiplication on the left 
by AA'j = a z A4ja z and subsequent diagonalization. The 
eigcnstates are then linear combinations of the type 



C\u x f\(r) , u± = J— 



(5) 
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where s x = sign(x), a = y/ j 2 — g 2 and f\(r) solves 

dlfx{r) + [e 2 + 2ge/r - a(a - A)/r 2 ] / A (r) = 0. (6) 

Introducing p = \e\r, the above becomes the familiar ra- 
dial equation for the 3D Coulomb problem [9], and the 
presence of e 2 (rather than e) entails the absence of bound 
solutions in the Dirac problem. It is important to note 
that when g is above g c = 1/2, the parameter a in eq. (6) 
becomes imaginary for some angular momentum chan- 
nels. The nature of the solutions is then radically differ- 
ent, which, as will be seen, has dramatic consequences. 
We address the two regimes separately. 

When g < g c , eq. (6) can be solved in terms of Coulomb 
wave functions [10, 11]: Fl(t],p), Gt{q,p). In fact, let- 
ting g = s e g, it is straightforward to show that the ap- 
propriate linear combination in (5) that solves (4) is 

(Pj(r)/J\fj = u + F a -!(-g, p) + s ge u_F a (-g, p), (7) 

where only the regular solution at the origin has been 
included. Since F a (—g,p) are the regular scattering so- 
lutions of the 3D Coulomb problem, they include the 
well-known logarithmic phase shift in the asymptotic ex- 
pansion [9]: 



F«(-9,p) ~ sin(p + glog(2p) +ti a (g) 



(8) 



where $ a (9) = — a§+arg[r(l-r-a— The logarithmic 
phase shift also carries to our case, for (7) can always be 
written asymptotically as 



p + glog{2p) + arg(u + e M ^ 1 



(9) 



The normalization, A/), is determined by imposing or- 
thogonality on the energy scale, J ipi(e,ryipj(e',r)dr = 
5 l0 8{e - e'), leading to A/T 2 = 2iT 2 a 2 /j 2 . With this 
choice, one conveniently recovers the free DOS per unit 
area and cone when g = 0. To see this, one notes 
that the LDOS, N{e,r) = T, E \^E(r)\ 2 S(e - E), 



is 



given by N(e,r) = ESL-oo n j( e > r )- Using rij(e,r) 



r 1 |< ( 9^ 1 (r)| 2 + r 1 |< y 9^(r)| 2 , the contribution from each 
angular momentum channel is simply: 

n 3 {e,r) = (N]/r) [F^ + F 2 + 2gF a F a ^/\j\] . (10) 

In the limit g — > the Coulomb wave functions reduce to 
Bessel functions [10], and one obtains N(e,r) = \e\/2n. 

Of the several aspects encoded in (10), two are im- 
mediate: particle-hole symmetry is lost, and the LDOS 
becomes singular as E — > 0. This last point follows from 
the fact that, in this limit, N(e,r) oc \e\ 2a and a < 1/2 
for \j\ = 1/2; the asymmetry stems from the dependence 
of g on the sign of the energy. 

It is most instructive to compare the results derived 
within the Dirac approximation (2), with the results on 
the lattice that one obtains using the full Hamiltonian in 




FIG. 1: (color online) (a) Comparison of the LDOS (solid) 
with the numerical results in the lattice (dashed), calculated 
at different distances from the impurity and g = 1/6. The 
DOS for g = is also included for comparison (dot-dashed). 
For clarity, curves at different r have been vertically displaced, 
(b) LDOS at the site closest to the impurity in the lattice for 
g = 1/3. (c) Quantization condition (11) with j — 1/2, for 
g = 0.1 and 0.4. 



eq. (1). In Fig. 1(a) one can observe that, at low ener- 
gies {E < 0.5i), the result of (10) reproduces the LDOS 
on the lattice even at distances of the order of the lat- 
tice parameter (the two cases are barely distinguishable 
for most of the plotted range). Moreover, the attractive 
Coulomb potential brings locally a reduction of spectral 
weight in the lower band, the opposite happening to the 
upper band. The effect is strongest near the impurity and 
evolves towards the bulk behavior at larger distances. 

This behavior of the spectrum near the Dirac point 
can be understood from an investigation of the quan- 
tized energies when the system is restricted to a region 
of finite radius, R. A convenient way to confine 2D Dirac 
fermions is to introduce an infinite mass at the boundary 
[12], which translates into the Boundary Condition (BC) 
(Pj(R) = tpf (R) . From (7) this expands into 

Qj = F a ^(-g,R\e\) - s jSt F a {^g,R\e\) = 0. (11) 

This equation, whose roots determine the quantized en- 
ergy levels, always has the trivial solution, e = 0. As one 
can easily verify [see Fig. 1(c) and also the asymptotic 
phase shift in (9)], the other nodes are simply shifted to 
lower energies with increasing g — just as expected un- 
der an attractive potential — but they never cross e = 0. 
This means that no states are phase-shifted across the 
Dirac point, but they rather heap up close to it when 
e > 0, and conversely for e < 0. Hence, even though 
the gap in graphene is zero, no states will cross it while 
g < g Cl much like in a conventional semiconductor. In the 
spirit of Friedel's argument [13], this effect has profound 
consequences for the induced charge, which will be dis- 
cussed later. Finally, we remark that, although the con- 



tinuum approximation does not support bound solutions 
(unless a cutoff is introduced), they appear naturally in 
the lattice. This is seen in Fig. 1(b), where a bound state 
barely detached from the band is signaled by the sharp 
peak at the lower band edge. Furthermore, notice how 
the van Hove singularities are strongly renormalized by 
the presence of the impurity. 

For g > g c , a can become purely imaginary, and we 
introduce [3 = — ia = \J g 2 — j 2 for those j's such that 
\g\ > \j\. In general, linearly independent solutions of 
(6) are [14] 

A = +l: F Q _i(-ff,p) andF_ Q (-g,p), (12) 
A = -l: F a (-g,p) and F_ a _ l {-~g,p). (13) 

When a€l the ones with negative index are divergent 
at the origin and thus only the first were kept in (7). 
But when a S iR, the solutions are well behaved at the 
origin (albeit oscillatory), and two linearly independent 
solutions emerge. One is analogous to (7): 

Cpifiir) = u+Fip-i(-g, p) + s jge u-Fip{-g, p), (14) 

apart from a normalization factor, and where now 



u± 




(15) 



The other solution is simply <^_^(r). The general so- 
lution is therefore of the type tfij(r) = C\ f>ip{r) + 
C2 fi-if)(r) , where C1.2 arc to be set by the BC at 
short distances. Since we seek the effective low energy 
description of a problem defined in a lattice, a natu- 
ral BC is to have an infinite mass at some short cut- 
off distance ao ~ a. This has the effect of forbidding 
the penetration of electrons to distances shorter than ao 
[12], thus reflecting the physical situation, while, at the 
same time, naturally curing the divergence in the po- 
tential at the origin. This translates again into a BC 
<pf{&o) = <fif (°o)j and given that Ci,2 can always be 
chosen so that C1/C2 = exp[2i<5j(e)], one then obtains 
the phase Sj(e): 



F. 



i/3— 1 



i/3 



i/3-1 



SejF lf3 



(16) 



We can follow the same procedure as before to normal- 
ize the states in the energy scale, and then extract the 
contribution of the overcritical j's to the LDOS: 



1 



g I J (p) + s tJ Re[e^g I /(p) 



n,j(e,r) 

2ir r (qj{oo) + s ej Re[e i2 ^f?f (00)] 
where, for readability, we defined 

q) = \F l3 \ 2 + \F lP ^\ 2 + ?M Re^F.^.x], 



(17) 



g 1 / = 2F i pF i0 _ 1 + l 4(F t 
9 
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FIG. 2: (color online) (a-c) The LDOS in the lattice (dashed, 
recursion method) is compared with the first contribution in 
(18) (solid) for different distances from the impurity, (d) The 
second contribution in (18) (solid) and the free, linear, DOS 
for reference, (e) First contribution in (18); the inset is a 
magnification for e ~ 0. In all panels g = 4/3. 



and (• • • ) r stands for the constant term as r — > 00. 
Eqs. (10) and (17) determine the LDOS for any coupling 
strength, g, which can be summarized as 



N(e,r) 



\o\<\g\ 



nj (e,r). (18) 



j\>\g\ 



The presence of the first term in eq. (18) brings a pro- 
found rearrangement of the spectrum close to the impu- 
rity, with much more striking consequences than in the 
weak coupling regime. In Fig. 2(a-c) we plot the LDOS 
obtained numerically in the lattice with g = 4/3, together 
with the first contribution in (18). For such g it comes 
only from ri±i/2(e, r), and we used ao = 0.55a to impose 
the BC. It is clear that the analytical result captures 
quite accurately the behavior of the LDOS in the lattice. 
Most importantly, both results exhibit 3 marked reso- 
nances in the negative (hole-like) energy region, which 
decay away from the impurity. Indeed, their amplitude 
is such that they dominate the profile of the LDOS at 
low energies. Increasing g will cause the resonances to 
migrate downwards in energy, and their number to in- 
crease. This is rather peculiar and has to do with the 
fact that, in reality, the Dirac point is an accumulation 
point of infinitely many resonances [inset in Fig. 2(e)]. 
One can appreciate the origin of this from the fact that 
F L -i(f},p~ 0) - p L . Since L G C in eqs. (14) and (16), 
it implies that the wave functions oscillate with logarith- 
mically diverging frequency as e — > 0. This situation is 
akin to the fall of a particle to the center [9] , and the ef- 
fect carries to the LDOS with the consequences shown in 
Figs. 2(a-c). In panel (d), we present the remainder con- 
tribution (second term) to the total LDOS in eq. (18). It 
is evident that in the region e < 0, dominated by the res- 
onances, this contribution is highly suppressed, whereas, 
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FIG. 3: (color online) (a) Induced charge numerically ob- 
tained from exact diagonalization on a lattice with 124 2 sites 
(g < g c ). (b) Idem (g > g c ). (c) Evolution of the numerical 
spectrum with g. (d) Analytical 8n(r) obtained using (19) 
and regularization. 

for positive energies, the LDOS exhibits an oscillating 
behavior around the bulk limit. 

As is customary, it is also of interest here to understand 
how the electronic density readjusts itself in the presence 
of this charged impurity. Even though interactions are 
not included (and thus there is no real screening), one 
can obtain important insights from the non-interacting 
problem in the spirit of Friedel [13] . The charge density, 
n(r), is straightforwardly obtained from integration of 
(18) in energy from an energy cut-off, —A, up to Ep = 0. 
Since it involves integrals of F a (—g, p), this can be done 
exactly. For example, when g < g c , one has (for each j 
channel) 

n 3 -(r) - [ P n 3 (e,r)/r ~ \j\F a F a ^/(7T 2 r 2 )\ e= _ A . (19) 

Although the above needs only to be evaluated at e = 
—A, it is not free from difficulties yet, for there is an 
infinite sum over j to be performed. Expanding (19) 
asymptotically, the induced charge per channel reads 

Sn J (r)=n ] (r)-n°Ar)^(l/r) [A-g/r-A + O(r- 2 )] , 

(20) 

where the remainder is oscillating with frequency A, and 
convergent with j. Clearly, if A = A the sum over j 



diverges. We regularize this by locally changing the cut- 
off: A = A + g/r, whereupon the leading contribution 
is ~ r -3 , and oscillating with frequency A [Fig. 3(d)]. 
Nonetheless, despite accidentally reproducing the lattice 
behavior, the oscillation itself is tied to the cutoff proce- 
dure. We point out that any charge oscillation decaying 
faster than 1/r 2 on the lattice appears, in the continuum 
theory, as a Dirac delta function at the origin, in agree- 
ment with perturbative studies of this problem [15], but 
differs from the self-consistent calculation in ref. [8]. In- 
terestingly, the behavior of the induced charge in the lat- 
tice is indeed ~ r -3 and oscillating, as seen in Fig. 3(a), 
wherein exact numerical results in the lattice are plotted. 
An analogous analytical procedure can be undertaken for 
g > g Cl leading to an induced charge decaying as ~ r~ 2 . 
Fig. 3(b) shows that this agrees with the numerical data 
in the lattice, where 5n(r) ~ ?'~ 2 , and non-oscillating. 
One thus concludes that the induced charge behaves quite 
differently below and above g c , as had been hinted before 
on account of the peculiar behavior of the phase shifts be- 
low g c . This last point can be confirmed by inspecting 
the behavior of the numerical energy levels as a func- 
tion of g shown in Fig. 3(c), being evident the difference 
between the two regimes. 

We have studied the problem of a Coulomb charge in 
graphene via exact numerical methods on the lattice and 
the Dirac Hamiltonian. We calculated the LDOS and 
local charge as a function of energy and distance from 
the impurity, having found that the Dirac equation pro- 
vides a qualitative description of the problem at low en- 
ergies. We found new features in the lattice description 
that are beyond the Dirac equation: bound states and 
strong renormalization of the van Hove singularities. We 
have also shown the existence of a critical coupling g c 
separating the weak and strong coupling regimes, with 
radical differences in the features of the LDOS. These 
results can be tested experimentally through STS mea- 
surements. 
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